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Four  technical  reports  covering  the  accomplishments  in  this  project  are  at¬ 
tached.  Below  we  provide  an  executive  summary  of  the  work  performed. 


Objectives 

•  Perform  a  Large  Eddy  Simulation  of  a  flow  around  an  airfoil  at  maxi- 
miun  lift  using  both  structured  and  imstructured  grids. 

•  Perform  efficient  parallel  implementation  of  the  codes. 

The  results  of  the  first  stage  of  the  investigations,  using  both  a  finite- 
difierence  method  on  a  structured  mesh  and  a  finite-element  method  on  an 
unstructured  mesh,  differed  considerably  from  each  other  and  from  the  avail¬ 
able  experimental  data.  Differences  were  found  with  respect  to  the  occurence 
of  transition  near  the  suction  peak  and  with  respect  to  the  amount  of  backflow 
near  the  trailing  edge.  To  identify  the  cause  of  the  mismatch,  we  addressed 
some  important  issues: 

•  The  sensitivity  of  the  overall  flow  to  details  in  the  nose  section,  in 
particular  to  the  transition  strip. 

•  The  effect  of  windtunnel  walls. 

•  The  spanwise  domain  size. 

To  answer  these  questions,  we  selected  the  finite-element  code  due  to  its 
ability  to  easily  handle  complex  geometry  and  local  refinement  of  the  mesh. 
At  the  same  time,  we  started  the  porting  of  the  structured  code  on  parallel 
machines. 


Accomplishment  and  future  plans 

Unstructured  code 

This  code  solves  the  compressible  Navier-Stokes  equations  on  an  unstructured 
grid  via  a  Galerkin/least-squares  stabilized  finite  element  method.  The  code 
has  been  ported  to  the  IBM  SP2.  The  code  now  uses  MPI,  a  communication 
standard  that  is  more  widely  used  than  that  of  the  original  code,  which  should 


facilitate  the  port  to  other  platforms.  As  indicated  in  the  attached  technical 
reports,  the  calculations  were  in  disagreement  with  the  experimental  data. 
Several  possible  causes  for  the  discrepancies  were  investigated. 

1.  Effect  of  wind  tunnel  walls  and  transition  strip 

In  the  experiment,  a  strip  of  tape  with  serrations  cut  into  the  edge  on  the 
upstream  side  was  used.  The  serrated  tape  has  been  modeled  in  a  coarse 
fashion  by  our  current  simulation  as  can  be  seen  in  Appendix  4,  Fig.  1.  The 
tape  is  effectively  a  forward  facing  step  (with  serrations)  of  height  <^99/4, 
followed  by  a  backward  facing  step. 

The  blockage  effect  of  the  wind  timnel  walls  has  also  been  included  in 
the  recent  calculations.  Note  that  the  boimdary  layers  on  the  walls  are  not 
simulated;  rather,  sUp  boundary  conditions  are  appUed  on  the  wind  tunnel 
walls  as  can  be  seen  in  Appendix  4,  Fig.  2. 

These  two  effects  were  studied  separately  for  a  short  period  of  time  (not 
sufficient  for  converged  statistics  in  the  trailing  edge  region),  and  agreement 
with  experiments  was  seen  to  improve  in  both  cases.  The  effect  of  the  walls 
was  somewhat  greater  than  that  of  the  transition  strip.  The  results  of  this 
simulation  can  be  seen  in  Appendix  4,  Fig.  3,  where  the  velocity  profiles  of  the 
new  simulation  (with  wind  timnel  walls  and  transition  strip)  are  compared  to 
the  original  simulation.  Note  the  large  increase  in  the  degree  of  separation. 
This  study  is  being  continued  by  Prof.  Jansen  at  Renssalear  Polytechnic 
Institute  under  a  separate  AFOSR  contract. 

2.  More  spanwise  domain  studies 

Based  on  the  findings  in  the  previous  section,  more  attention  was  given  to 
the  spanwise  domain  effects.  Since  the  resolution  changed  at  the  same  time 
as  the  expansion  of  the  domain,  it  is  difficult  to  isolate  the  two  effects.  For 
this  reason  the  first  study  maintained  5%  chord  domain  width  and  improved 
the  resolution  in  the  second  half  of  the  airfoil  where  the  solution  has  shown 
some  change.  These  calculations  are  being  carried  out  by  Prof.  Jansen. 

3.  Higher  order  methods 

Given  the  number  of  points  that  are  required  to  obtain  a  grid-independent 
solution,  it  seems  clear  that  higher  order  methods  should  be  explored.  This  is 
straightforward,  but  non-trivial,  to  do  with  the  finite  element  method.  There 


is  another  benefit  to  higher  order  methods  besides  the  obvious  one  of  higher 
accuracy.  The  higher  order  methods  will  have  a  more  complete  representation 
of  the  residual  error  of  the  discrete  approximation  and,  therefore,  the  scheme 
will  be  less  dissipative. 

A  fourth-order  finite  difference  code  with  suitable  conservation  properties 
has  also  been  developed.  This  code  has  been  tested  in  turbulent  channel  flow. 

Structured  code 

The  code  is  based  on  a  second-order  finite  difference  method  which  solves,  on 
a  staggered  mesh,  the  three-dimensional  incompressible  Navier-Stokes  equa¬ 
tions  in  primitive  variables  (velocity  and  pressure)  and  in  generalized  coordi¬ 
nates  on  a  spanwise  periodic  domain.  The  solution  is  advanced  in  time  by  a 
spTni-implir.it  third-order  Runge-Kutta  scheme.  The  SGS  model  implemented 
is  a  version  of  the  dynamic  model  suitable  for  application  in  generalized  co¬ 
ordinates,  with  the  coefficient  averaged  in  the  spanwise  direction.  The  use 
of  a  C-mesh  allows  a  better  resolution  in  the  wake  region  and  simplifies  the 
imposition  of  outflow  boundary  conditions. 

The  serial  code  was  optimized  on  a  Cray  C90.  Different  strategies  should 
be  used  on  different  processors  found  in  parallel  computers,  where  perfor¬ 
mance  is  heavily  affected  by  the  presence  of  memory  hierarchies.  The  array 
layout  has  been  rearranged  and  the  code  has  been  optimized  on  a  single  SP2 
processor. 

The  parallel  formulation  employs  spatial  decomposition  of  the  overall 
grid  into  sub  blocks  that  are  assigned  to  separate  processors  and  use  MPI  as 
the  message-passing  programming  interface.  A  2-D  domain  decomposition 
in  the  plane  of  generalized  coordinates  is  used  in  order  to  distribute  data 
across  concurrent  processors.  The  2-D  decomposition  will  allow  the  use  of  a 
large  number  of  processors  (64,128)  with  a  good  ratio  between  surfaces  and 
volumes  of  the  domains.  A  flexible  decomposition  has  been  implemented 
in  which  each  domain  could  have  a  different  size  (the  2-D  decomposition 
is  the  product  of  two  general  1-D  decompositions)  and  different  boundary 
conditions. 

The  solution  at  points  shared  by  neighboring  processors  is  updated  be¬ 
tween  each  sub  step  by  means  of  a  message  exchange.  Due  to  the  staggered 
grid  and  generalized  coordinates,  each  block  needs  two  layers  of  “ghost”  cells. 
Some  important  issues  such  as  the  solution  of  the  Poisson  equation  necessary 
to  project  the  velocity  field  in  a  divergence-free  space  have  been  addressed. 


A  parallel  multigrid  solver  has  been  implemented  using  MPI  as  the  commu¬ 
nication  hbrary.  This  solver  utilizes  a  parallel  tridiagonal  solver,  and  it  is 
capable  of  full-coarsening  and  semi-coarsening. 

The  parallel  code  will  permit  the  use  of  fine  grids  that  will  avoid  unphys¬ 
ical  oscillation  observed  when  the  mesh  is  too  coarse. 
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APPENDIX  1.  Preliminary  large-eddy 
simulations  of  flow  around  a  NACA 
4412  airfoil  using  unstructured  grids 

By  Kenneth  Jansen 


1.  Motivation  and  objectives 

Large-eddy  simulation  (LES)  has  matured  to  the  point  where  application  to  com¬ 
plex  flows  is  desirable.  The  extension  to  higher  Reynolds  numbers  leads  to  an  im¬ 
practical  number  of  grid  points  with  existing  structured-grid  methods.  Furthermore, 
most  real  world  flows  are  rather  difficult  to  represent  geometrically  with  structmed 
grids.  Unstructured-grid  methods  offer  a  release  from  both  of  these  constraints. 
However,  just  as  it  took  many  years  for  structured-grid  methods  to  be  well  un¬ 
derstood  and  reliable  tools  for  LES,  vmstructured-grid  methods  must  be  carefully 
studied  before  we  csm  expect  them  to  attain  their  full  potential. 

In  the  past  two  years,  important  building  blocks  have  been  put  into  place  mak¬ 
ing  possible  a  careful  study  of  LES  on  unstructured  grids.  The  first  building  block 
was  an  efficient  mesh  generator  which  allowed  the  placement  of  points  according  to 
smooth  variation  of  physical  length  scales.  This  variation  of  length  scales  is  in  all 
three  directions  independently,  which  allows  a  large  reduction  in  points  when  com¬ 
pared  to  structured-grid  methods,  which  can  only  vary  length  scales  in  one  direction 
at  a  time.  The  second  building  block  was  the  development  of  a  dynamic  model  ap¬ 
propriate  for  unstructured  grids.  The  principle  obstacle  was  the  development  of 
an  imstructiured-grid  filtering  operator.  New  filtering  operators  were  developed  in 
Jansen  (1994).  In  the  past  year,  some  of  these  filters  have  been  implemented  into 
a  highly  parallelized  finite  element  code  based  on  the  Galerkin/least-squares  finite 
element  method  (see  Jansen  tt  al.  (1993)  and  Johan  et  al,  (1992)). 

We  have  chosen  the  NACA  4412  airfoil  at  maximum  lift  as  the  first  simulation 
for  a  variety  of  reasons.  First,  it  is  a  problem  of  significant  interest  since  it  would 
be  the  first  LES  of  an  aircraft  component.  Second,  this  flow  has  been  the  subject 
of  three  experimental  studies  (Coles  and  Wadcock  (1979),  Hasting  and  Williams 
(1987),  and  Wadcock  (1987)).  The  first  study  found  the  maximum  lift  angle  to  be 
13.87°.  The  later  studies  found  the  angle  to  be  12°.  Wadcock  reports  in  the  later 
study  that  the  early  data  agree  very  well  with  his  new  data  at  12°,  suggesting  that 
the  early  experiment  suffered  from  a  non-parallel  mean  flow  in  the  Caltech  wind 
tunnel.  It  should  be  pointed  out  that  the  Reynolds-averaged  simulations  are  usually 
run  at  13.87°  and  do  not  agree  with  the  data  when  rim  at  12°  as  will  be  shown  later 
in  this  study.  It  is  hoped  that  LES  can  clarify  this  controversy.  The  third  reason  for 
considering  this  flow  is  the  variety  of  flow  features  which  provide  an  important  test  of 
the  dynamic  model.  Starting  from  the  nose  where  the  flow  stagnates,  thin  laminar 
boundary  layers  are  formed  in  a  very  favorable  pressure  gradient.  This  pressure 
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gradient  soon  turns  adverse,  driving  the  flow  toward  a  leading  edge  separation. 
Only  the  onset  of  turbulence  can  cause  the  flow  to  remzdn  attached  or  to  reattach 
if  it  did  separate.  The  persistent  adverse  pressure  gradient  eventually  drives  the 
turbulent  flow  to  separate  in  the  last  20  percent  of  chord.  The  separation  bubble  is 
closed  near  the  trailing  edge  as  the  reteirded  upper  surface  boundary  layer  intereicts 
with  the  very  thin  lower  surface  boundary  layer.  The  large  difference  in  boundary 
layers  creates  a  challenging  wake  to  simulate.  Only  the  dynamic  model  can  be 
expected  to  perform  satisfactorily  in  this  variety  of  situations:  from  the  laminar 
regions  where  it  must  not  modify  the  flow  at  all  to  the  turbulent  boimdary  layers 
and  wake  where  it  must  represent  a  wide  variety  of  subgrid-scale  structures. 

The  flow  conflguration  we  have  chosen  is  that  of  Wadcock  (1987)  at  Reynolds 
number  based  on  chord  Rcc  =  Uooc/u  =  1.64  x  10®,  Madi  number  M  =  0.2,  and 
12®  angle  of  attack. 

2.  Accomplishments 

2.1  Dynamic  model  implemented  and  tested 

The  only  obstacle  to  implementing  a  dynamic  model  on  unstructured  grids  is 
extension  of  the  filtering  operator.  Four  filtering  operators  were  proposed  in  Jansen 
(1994).  Two  of  these  models  were  implemented  and  compared  using  a  simple  ana¬ 
lytic  velocity  field  for  which  the  filtered  values  can  be  determined  exactly.  Prom  this 
test,  the  generalized  top-hat  was  found  to  be  the  most  accurate,  and  all  subsequent 
calculations  have  been  carried  out  using  this  filter. 

2.2  Simulations 

A  series  of  simulations  has  been  performed  in  the  last  year  to  develop  experi¬ 
ence  with  this  new  approach.  The  first  simulation  was  intentionally  very  coarse 
as  we  hoped  to  improve  the  mesh  selectively  and  develop  an  understanding  of  the 
sensitivity  of  the  solution  to  the  grid  improvements. 

2.2.1  First  simulation 

The  first  simulation  was  performed  on  a  very  coarse  mesh.  The  near-wall  grid 
did  not  attempt  to  resolve  the  near-wall  layer  accurately  in  the  first  20  percent  of 
chord  and  only  marginally  resolved  the  remaining  flow  (A^  =  300,  Aj  =  80)  at 
the  wall.  The  grid  was  coarsened  in  the  streamwise  and  spanwise  directions  coming 
off  the  wall  as  suggested  by  Chapman  (1979).  The  resulting  coefficient  of  pressiure 
distribution  was  reasonably  well  predicted  on  this  mesh  (see  Fig.  1),  but  the  velocity 
profiles  showed  poor  agreement  with  the  experiment. 

2.2.2  Improvement  of  outer  layer 

Careful  scrutiny  of  the  mesh  revealed  that  the  strategy  of  coarsening  in  the 
streamwise  direction  coming  off  the  wall  was  inappropriate  at  this  Reynolds  number. 
The  inner-layer  spacing.  A,  scales  on  wall  xmits. 


APPENDIX  1.  LES  with  unstructured  grids 
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Figure  1.  Coefl5.cient  of  pressure  along  the  airfoil  surface.  LES  —  ,  Hastings 

and  WilUams  o  ,  Wadcock  (1987)  a  . 

Here,  v,  is  the  kinematic  viscosity,  c  is  the  chord  length,  and  Mr  is  the  friction 
velocity  defined  to  be  the  square  root  of  the  coefficient  of  friction,  C/,  over  two. 
The  outer-layer  spacing  scales  on  the  boundary  layer  thickness,  ^gg.  It  is  reasonable 
to  expect  the  large  eddies  in  the  outer  part  of  the  boundary  layer  to  be  of  order  ^gg, 
and  therefore  the  outer-layer  spacing,  in  all  directions,  should  never  exceed 


By  nsingr  Wadcock’s  experimental  data  for  the  Cf  and  ^gg,  one  can  compare  these 
two  resolution  restrictions  as  is  done  in  Fig.  2.  This  figure  contains  three  curves. 
The  solid  curve  describes  the  variation  of  a  200  wall-unit  spacing  (which  can  be 
associated  with  the  streamwise  spacing  near  the  wall)  over  the  upper  surface  where 
the  boundary  layer  is  attached.  The  dashed  curve  describes  the  same  variation 
of  50  wall  Tinits  (which  can  be  associated  with  spanwise  spacing  near  the  wall). 
The  chain  dash  curve  is  the  outer-layer  spacing  as  described  above.  Several  points 
can  be  made  in  this  figure.  First,  all  three  curves  change  by  over  an  order  of 
magnitude  from  the  tip  to  the  tail  region.  This  illustrates  how  an  unstructured  grid 
saves  points  by  matching  resolution  to  the  local  changes  in  the  length  scales  in  the 
streamwise  direction.  For  example,  a  structured  grid  would  be  forced  to  carry  the 
fine  spanwise  resolution  required  near  the  nose  through  the  entire  domain.  Second, 
when  comparing  the  near-wall  spanwise  resolution  to  the  outer-layer  resolution, 
it  is  clear  that  coarsening  the  spanwise  resolution  as  the  distance  from  the  wall 
increases  is  justified.  The  final  point,  apparent  from  this  figure,  is  that  coarsening 
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Figure  2.  Comparison  of  length  scales  over  the  airfoil  surface.  200  wall  xmits 

(streamwise  near- wall  spacing) - ,  ^99  /5  (outer-layer  spacing) - ,  50  wall 

maits  (spanwise  near-wall  spacing) - . 

of  the  streamwise  resolution  in  the  outer  layer  is  not  justified.  In  fact,  over  much 
of  the  airfoil  surface  the  outer-layer  grid  resolution  is  more  restrictive  than  the 
inner-layer  resolution.  The  choices  of  200  wall  units  and  5  points  per  boamdaxy 
layer  thickness  Me  somewhat  arbitrary,  but  they  are  believed  to  be  comparable  in 
their  degree  of  coarseness.  It  is  interesting  to  observe  that  the  crossover  between 
these  two  curves  corresponds  to  Rcr  =  (ur^as/*')  =  1000.  Therefore,  when  above 
1000,  the  inner-layer  resolution  is  the  most  restrictive.  Otherwise,  the  outer-layer 
resolution  is  the  most  restrictive.  Only  at  higher  Reynolds  numbers  will  coarsening 
in  the  streamwise  direction  be  justified. 

Considering  the  above  discussion,  a  new  mesh  was  made  where  the  coarsening 
of  the  streamwise  spacing  was  delayed  until  outside  of  the  boundary  layer.  This 
resulted  in  a  mesh  with  nearly  twice  as  many  points  as  the  previous  simulation.  It 
also  resulted  in  a  rather  dramatic  change  in  the  early  boundary  layer  structme.  It 
seems  that  the  improved  resolution  of  the  outer  layer  allowed  a  better  resolution  of 
the  leading  edge  separation.  The  new  simulation  led  to  a  train  of  spanwise  coherent 
vortices.  These  vortices  broke  down  into  turbulence  at  about  10  percent  of  chord. 

The  persistence  of  the  spanwise  coherent  vortices  was  not  in  line  with  the  exper¬ 
iments  which  were  all  tripped.  Some  evidence  as  to  the  importance  of  the  tripping 
can  be  seen  in  Fig.  3  where  we  compare  the  surface  coefl&cient  of  pressure  distribu¬ 
tion  from  the  free  transition  simulation  to  two  experimental  data  sets  from  Hastings 
and  Williams  (1987).  The  square  data  set  was  taken  without  a  transition  strip  while 
the  circle  data  set  was  taken  with  a  transition  strip.  Ovur  simulation  shows  rather 
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good  agreement  with  the  free  transition.  Unfortunately,  all  velocity  and  Rejmolds 
stress  data  were  taJcen  with  the  transition  strip  in  place  and  agreement  with  these 
quantities  is  substantially  worse  them  with  Cp. 


Figure  3.  Coeflficient  of  pressure  along  the  airfoil  surface.  LES - ,  Hastings 

and  Williams  without  trip  □  ,  Hastings  and  Williams  tripped  o  . 

2.2.S  Grid  refinement  study  of  the  nose 

The  dramatic  change  of  the  flow  with  changed  resolution  indicated  a  need  for 
further  refinement  in  the  nose  re^on.  At  the  same  time  we  also  hoped  to  model  the 
transition  through  a  steady  blowing  pattern  as  shown  in  Fig.  4.  A  shape  that  could 
be  easily  resolved  was  chosen.  Therefore,  we  could  be  certain  that  any  sensitivity 
to  grid  refinement  would  be  associated  with  the  tmrbulence  structures  responding 
to  the  blowing  and  not  the  resolution  of  the  blowing  itself. 

A  new  mesh  was  generated  where  the  streamwise  and  spanwise  resolution  were 
improved  by  a  factor  of  two  everywhere  on  the  upper  surface.  The  normal  spacing 
was  improved  at  the  wall  by  a  factor  of  two  as  well,  but  this  did  not  lead  to  a 
doubling  of  points  in  this  direction  due  to  the  stretching.  The  spanwise  domain 
was  cut  in  half  (from  0.05c  to  0.025c)  for  this  simulation.  Therefore,  the  number  of 
points  approximately  doubled  rather  than  a  quadrupling. 

There  was  again  a  rather  dramatic  change  in  the  solution  and  so  another  mesh 
was  generated.  This  mesh  again  improved  the  streamwise  and  spanwise  resolution 
by  a  factor  of  two,  although,  this  time,  only  in  the  first  5  percent  of  chord.  The  three 
surface  meshes  of  the  first  10  percent  of  chord  are  shown  in  Fig.  5.  The  velocity 
profiles  in  the  first  5  percent  of  chord  are  shown  in  Fig.  6.  For  this  forcing  pattern, 
the  flow  is  nearly  spanwise-  and  streamwise-resolution  independent. 
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Figure  4.  Elevation  plot  of  the  steady  jet  normeil  to  the  airfoil  surface.  The 
actual  grid  is  shown  to  confirm  resolution. 


Figure  5.  Sxuface  meshes  near  the  leading  edge  (0.0  <  x/c  <  0.1).  Mesh  (b)  has 
been  refined  by  a  factor  of  2,  both  spanwise  and  streamwise,  from  mesh  (a).  The 
spanwise  domain  is  also  halved.  Mesh  (c)  has  been  refined  by  a  factor  of  2,  both 
spanwise  and  streamwise,  from  mesh  (b)  in  the  first  5  percent  of  chord. 
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Figure  6.  Profiles  of  tangential  velocity  component  at  various  positions  along 
the  airfoil  surface  (x/c  =  0.005,0.01,0.02,0.03,0.04,0.05).  Solutions  correspond  to 
grids  firom  previous  figure;  mesh  (a) - ,  mesh  (b)  ,  mesh  (c)  . 


Figure  7.  Boundary  layer  parameters  (lower  curves  and  data  are  momentum 
thi^Vripss  (0),  higher  cmrves  and  data  are  displacement  thickness  (6*))  along  the 

airfoil  surface.  Solutions  correspond  to  meshes  from  Fig.  5;  mesh  (a) - ,  mesh 

(b) - ,  mesh  (c) - ,  Wadcock  (1987)  (□  ). 
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2.3  More  accurate  transition 

While  it  was  useful  to  obtain  a  grid-independent  solution  at  the  forcing  prescribed, 
the  final  solution  does  not  agree  with  experiment,  as  can  be  seen  in  Fig.  7.  Here, 
the  momentum  and  displacement  thickness  of  the  grid-independent  calculations 
can  be  seen  to  be  substantially  greater  than  the  experiments  at  the  first  available 
datum  point.  The  discrepancy  seems  to  be  associated  with  the  laminar  separation 
at  1  percent  of  chord.  The  simulation  suggests  a  transition  in  the  free  shear  layer, 
followed  by  a  turbtilent  reattachment.  This  mode  of  transition  seems  to  give  the 
flow  a  large  jump  in  momentum  and  displacement  thickness.  The  experiment  did 
not  seem  take  this  route  to  transition.  For  this  reason  a  more  careful  study  of 
transition  is  currently  underway. 

Wadcock  used  a  strip  of  tape  with  serrations  cut  into  the  edge  on  the  upstream 
side.  The  serrated  tape  can  be  modeled  in  a  coarse  fashion  by  our  current  simulation 
as  can  be  seen  in  Fig.  8.  The  tape  is  effectively  a  forward  facing  step  (with  serrations) 
of  height  ^99/4,  followed  by  a  backward  facing  step.  Calculations  are  underway  with 
this  modification. 


Figure  8.  A  transition  strip  is  modeled  geometrically  by  applying  a  no-slip 
botmdary  condition  to  the  nodes  which  form  a  surface  of  height,  shape,  and  position 
equivalent  to  Wadcock’s  serrated  tape  which  was  applied  to  the  airfoil  surface. 

2.4  Reynolds-averaged  simulations 

Reynolds-averaged  Navier-Stokes  simulations  (RANS)  have  not  shown  good  agree¬ 
ment  with  the  experimental  data.  However,  given  the  cost  of  LES,  they  can  be  a 
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helpful  tool  for  suggesting  sensitivity  to  changes  of  basic  flow  parameters  since  they 
require  so  little  computational  time.  While  the  results  are  not  expected  to  be  quan¬ 
titatively  correct,  trends  can  at  least  be  suggested  by  RANS  and  later  confirmed 
by  LES. 

A  series  of  RANS  calculations  was  performed  to  chart  various  trends  in  this 
flow.  The  RANS  calculations  used  the  commonly  accepted  NASA  code  (INS2D)  of 
Rogers  (1991)  and  employed  a  A:  -  w  model  from  Menter  (1994).  First,  the  effect 
of  angle-of-attack  and  wind  ttmnel  walls  are  compared  in  Figures  9  and  10.  The 
boundary  condition  on  the  wind  tmrnel  walls  is  a  slip  condition.  This  accounts  for 
the  blockage  of  the  walls  without  requiring  resolution  of  the  boundary  layers  on 
them.  The  effects  are  compared  together  because  it  is  common  among  the  RANS 
modeling  community  to  adjust  the  angle-of-attack  of  free  air  calculations  to  account 
for  the  walls.  Figure  10  suggests  that  the  flattening  of  the  Cp  near  the  trailing  edge 
(which  is  associated  with  the  large  separation  there)  is  affected  strongly  by  angle- 
of-attack  and  only  weakly  by  the  wind  tunnel  walls.  The  13.87®  angle-of-attack 
cannot  be  justified  with  the  hope  of  accounting  for  the  effects  of  the  wind  tunnel 
walls  in  free  air  calculations. 

The  second  trend  studied  with  the  RANS  code  was  the  effect  of  transition  posi¬ 
tion.  When  the  RANS  code  was  rim  with  the  transition  point  fixed  at  the  position 
of  Wadcock’s  strip,  a  leading  edge  separation  developed  on  sufficiently  fine  meshes. 
Once  beyond  the  transition  point,  the  flow  reattaches.  This  provides  an  independent 
verification  of  the  results  observed  in  the  LES. 

3.  Future  plans 

S.l  Grid-independent  solution  of  flow  with  a  transition  strip 

The  calculation  using  the  transition  strip  described  above  will  be  continued  and 
checked  for  grid  dependence.  It  should  be  noted  that  grid  independence  can  only  be 
achieved  beyond  a  short  distance  downstream  of  the  transition  strip.  True  grid  inde¬ 
pendence  of  the  strip  and  transition  itself  is  probably  too  expensive  to  be  practical, 
even  with  an  unstructured  grid.  It  may  be  necessary  to  provide  small  disturbances 
upstream  of  the  strip  to  mimic  the  interaction  of  freestream  turbulence  with  the 
strip.  This  capability  has  been  implemented  and  tested  in  the  code  using  a  wall  jet 
with  spatial  and  temporal  variation. 

S.2  Inclusion  of  the  wind  tunnel  walls 

The  RANS  studies  indicated  a  moderate  effect  of  the  wind  tunnel  walls  on  the 
solution.  Future  simulations  will  be  done  with  a  slip  boundary  conditions  on  the 
wind  tunnel  walls.  Meshes  have  already  been  generated  for  this  purpose  as  can  be 
seen  in  Fig.  11. 

3.3  Higher  order  methods 

Given  the  number  of  points  that  are  required  to  get  a  grid-independent  solution, 
it  seems  clear  that  higher  order  methods  should  be  explored.  This  is  relatively  easy, 
but  non-trivial,  to  do  with  the  finite  element  method.  There  are  two  benefits  to 
higher  order  methods  besides  the  obvious  one  of  higher  accmacy.  First,  the  higher 


Figure  9.  Coefficient  of  pressure  along  the  airfoil  surface  from  RANS  simulations. 

12®  in  free  air - ,  12®  with  walls - ,  13.87®  in  free  air - ,  13.87®  with 

walls - ,  Hastings  and  WiUiams  (1987)  o  ,  Wadcock  (1987)  □  . 
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Figure  10.  Closeup  of  the  trailing-edge  region  (same  legend  as  above). 
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Figure  11.  The  crosssectional  plane  of  an  imstructured  mesh  which  accoimts  for 
the  inviscid  effects  of  the  wind  tunnel  walls. 

order  methods  will  have  a  more  complete  representation  of  the  residual  error  of  the 
discrete  approximation  and,  therefore,  the  scheme  will  be  less  dissipative.  Second, 
alternative  filters,  described  in  Jansen  (1994),  can  be  implemented  and  tested.  It 
is  difficult  to  predict  at  this  time  if  the  method  will  lose  computational  efficiency 
when  extended  to  higher  order. 

3.4  Expanded  spanwise  domain 

Once  we  are  satisfied  with  the  solution  in  the  region  of  the  nose,  we  will  have  to 
consider  carefully  the  effect  of  the  narrow  spanwise  domain  on  our  solution.  It  is 
likely  that  as  we  predict  a  larger  separation  at  the  trailing  edge,  the  effect  of  the 
narrow  domain  will  become  more  acute.  Strategies  are  being  developed  to  expand 
only  the  portion  of  the  domain  suffering  from  a  narrow  box.  If  these  strategies  work 
a  large  number  of  points  can  be  saved. 
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APPENDIX  2.  Large-eddy  simulation  of 
flow  around  an  airfoil  on  a  structured  mesh 

By  Hans- Jakob  Kaltenbach  AND  Haecheon  Choi 

1.  Motivation  and  objectives 

The  diversity  of  flow  characteristics  encountered  in  a  flow  over  an  airfoil  near 
mnYimuTn  lift  tajces  the  presently  available  statistical  turbulence  models.  This  work 
describes  our  flrst  attempt  to  apply  the  technique  of  large-eddy  simulation  to  a 
flow  of  aeronautical  interest.  The  challenge  for  this  simulation  comes  from  the  high 
Reynolds  number  of  the  flow  as  well  as  the  variety  of  flow  regimes  encountered, 
including  a  thin  laminar  boundary  layer  at  the  nose,  transition,  boimdary  layer 
growth  under  adverse  pressure  gradient,  incipient  separation  near  the  trailing  edge, 
and  merging  of  two  shear  layers  at  the  trailing  edge. 

The  flow  configuration  chosen  is  a  NACA  4412  airfoil  near  maximum  hft.  The 
corresponding  angle  of  attack  was  determined  independently  by  Wadcock  (1987) 
and  Hastings  &  Williams  (1984,  1987)  to  be  close  to  12°.  The  simulation  matches 
the  chord  Reynolds  number  Uooclv  =  1.64  x  10®  of  Wadcock’s  experiment. 

2.  Accomplishments 

2.1  Numerical  method  and  SGS  model 

The  numerical  method  for  solving  the  unsteady,  incompressible  Navier-Stokes 
equations  is  described  in  Choi  et  al.  (1993).  Second-order  spatial  central  differences 
on  a  staggered  mesh  are  combined  with  a  semi-implicit  time  integration  scheme. 
Formulation  of  the  problem  in  terms  of  contravariant  velocity  components,  weighted 
with  the  Jacobian,  in  conjimction  with  the  staggered  variable  configuration  leads  to 
discretized  equations  that  can  be  solved  with  the  classical  splitting  approach.  The 
resulting  pressure  Poisson  equation  is  solved  using  FFT  for  the  spanwise  (periodic) 
direction  and  iterative  methods  for  the  remaining  two-dimensional  problems.  The 
computational  cost  is  about  equally  distributed  between  computation  of  the  right- 
hand  side  and  solving  the  Poisson  equation  at  every  substep  of  a  third  order  Runge 
Kutta  time  integration. 

The  implementation  of  the  dynamic  subgrid-scale  model  (Germano  et  al.  1991) 
with  least-square  contraction  (Lilly  1992)  uses  the  spanwise  homogeneity  of  the  flow 
to  obtain  a  model  coefidcient  that  is  a  function  of  streamwise  and  wall-normal  coor¬ 
dinate  only.  We  found  that  the  dynamic  procedure  occasionally  renders  unrealistic 
negative  coefldcients  in  regions  where  the  flow  is  laminar  sudi  as  at  the  nose  or  in 
the  potential  flow  region.  In  these  regions,  the  negative  values  are  the  result  of  the 
dynamic  model  becoming  ill-conditioned  and  have  no  physical  sigmficance.  In  the 
present  simtdations  we  prevented  any  form  of  backscatter  by  constraining  the  model 
coefficient  to  be  always  positive. 
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x/c 

^99 /c 

Ax/c 

Ax'*" 

Az+ 

Ax/Az 

h99fl'z 

mm 

405 

137 

2.96 

mssM 

bB 

378 

118 

3.2 

wsm 

274 

86 

3.2 

0.6 

0.0050 

235 

49 

4.8 

0.6 

0.8 

■MW 

0.0088 

110 

13 

8.5 

1.2 

Table  1.  Spacing  along  upper  surface  of  airfoil.  The  last  two  columns  show  cell 
aspect  ratio  and  ratio  of  boundary  layer  thickness  to  domain  width  for  case  A. 


On  the  present  mesh,  the  CFL  limit  of  1.5  results  in  an  average  timestep  of 
2  X  10~^c/17oo.  About  80  CPU-seconds  on  a  Cray-C90  are  needed  to  advance  the 
solution  over  one  timestep  on  a  mesh  of  638  x  79  x  48  =  2.4  x  10®  cells.  Therefore, 
simulation  of  one  time  unit  cjUoo  requires  90  CPU-hours.  In  order  to  obtain  smooth 
statistics  the  results  have  to  be  averaged  over  several  time  units. 

2.2  Computational  domain  and  mesh  layout 

The  computational  domain  is  a  C-mesh  with  the  outer  boundary  about  three 
chord  lengths  away  from  the  surface.  At  the  outer  bovmdary  we  specify  the  freestream 
velocity  ?7oo.  As  a  consequence,  the  vertical  velocity  component  (in  a  coordinate 
system  aligned  with  the  chord  at  0®  angle  of  attack)  will  be  zero  at  the  outer  bound¬ 
ary.  Therefore,  the  chosen  configuration  resembles  more  the  flow  around  an  airfoil 
inside  of  a  wind  tunnel  with  parallel  walls  than  an  airfoil  in  free  air.  Jansen  (1995) 
has  shown  that,  even  with  the  walls  located  much  closer,  the  presence  of  wind  tun¬ 
nel  walls  mainly  affects  the  flow  in  the  nose  region  by  increasing  the  suction  peak. 
The  pressure  distribution  in  the  rear  part  and  the  size  of  the  backflow  zone,  how¬ 
ever,  are  only  weakly  dependent  on  whether  the  wind  tunnel  walls  are  included  or 
not.  A  no-sHp  condition  is  enforced  at  the  airfoil  surface,  and  we  use  a  convective 
(radiative)  boimdary  condition  at  the  outflow  plane. 

Results  from  two  simulations  will  be  presented.  The  two  cases  differ  only  with 
respect  to  the  spanwise  domain  width  which  is  0.05c  in  case  A  and  0.025c  in  case 
B.  The  spanwise  spacing  Aar  is  the  same  with  48  cells  in  case  A  and  24  cells  in 
case  B,  respectively.  Main  criterion  for  the  choice  of  the  spanwise  domain  size 
is  the  ratio  of  botmdary  layer  thickness  to  domain  width,  which  is  tabulated  in 
Table  1.  As  a  consequence  of  the  rapid  growth  of  the  boundary  layer  thickness  on 
the  suction  side,  this  ratio,  which  is  initially  sufficiently  small  to  capture  several 
structures  in  the  spanwise  direction,  exceeds  one  near  the  trailing  edge.  It  is  likely 
that  the  development  of  flow  structures  in  the  outer  part  of  the  boundary  layer  will 
be  affected  by  the  limited  domain  size.  Comparison  of  cases  A  and  B  gives  insight 
in  the  sensitivity  of  the  simulation  with  respect  to  this  parameter. 

The  design  of  an  adequate  mesh  involves  several  aspects.  The  most  energetic 
eddies  of  the  boundary  layer  have  to  be  resolved.  More  or  less  general  criteria 
have  been  developed  for  the  mesh  spacing  in  the  case  of  wall  boimded  shear  flows 
under  zero  pressure  gradient.  However,  these  criteria  depend  on  the  munerical 
method  employed  (Limd  et  al,  1995).  Cabot  (1994)  found  that  for  LES  of  turbvJent 
channel  flow  based  on  second-order  finite  differences  a  spacing  of  Ar'*"  =  60  and 
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(bottom)  and  xjc^  0.98  (top)  at  about  5%  of  the  local  botmdary  layer  height. 
Individual  curves  are  separated  by  a  vertical  offset  of  0.3  with  the  corresponding 
zero-lines  located  at  1.2, 1.5,  ...  3.9. 

=  15  —  20  is  needed  to  adequately  resolve  the  near  wall  structures. 

Little  is  known  about  the  minimum  spacing  requirements  for  boundary  layers 
which  are  close  to  separation.  The  mesh  size  in  terms  of  wall  imits  probably  becomes 
less  relevant  in  this  case.  About  half  of  the  640  streamwise  points  were  distributed 
over  the  upper  surface,  which  guaranteed  that  the  streamwise  spacing  was  between 
1/3  and  1/5  of  the  local  boimdary  layer  thickness  for  most  of  the  upper  surface,  see 
Table  1.  The  streamwise  spacing  varies  considerably  along  the  surface  due  to  the 
boundary  layer  growth.  Near  the  trailing  edge,  the  grid  was  refined  in  x  in  order 
to  resolve  the  merging  of  the  two  shear  layers.  No  attempt  was  made  to  resolve  the 
turbulence  on  the  lower  side  of  the  airfoil.  Spacings  in  terms  of  wall  units  based 
on  the  local  skin  friction  as  given  in  Wadcock’s  experiment  are  ^ven  in  Table  1.  It 
is  evident  that  the  spacing  in  the  present  simulation  is  considerably  coarser  than 
what  has  been  found  to  be  necessary  for  channel  flow  simulations.  However,  as 
the  boimdary  layer  develops  along  the  surface,  the  resolution  criteria  become  less 
restrictive  so  that  the  flow  in  the  rear  part  is  much  better  resolved  than  in  the  front 
section. 

In  the  wall-normal  direction  we  used  a  hyperbolic  mesh  generator  (Chan,  1993) 
to  distribute  79  layers  of  cells.  The  first  line  away  from  the  wall  was  at  about 
y+  =  1,  and  over  most  of  the  surface  there  were  between  20  and  30  points  inside 
the  boundary  layer. 

2.S  Difficulties  arising  from  the  high  Reynolds  number 

Centered  difference  schemes  suffer  from  the  emergence  of  grid-to-grid  oscillations 
(2-A-waves,  wiggles)  when  used  for  high  Reynolds  number  simulations.  Usually, 
the  viscosity  provided  by  the  subgrid-scale  model  is  sufficient  to  dampen  these  grid- 
to-grid  oscillations.  Several  somrces  for  2-A-waves  have  been  identified  in  the  past 
(Gresho,  1981).  They  include  high  cell  Peclet  numbers  in  conjunction  with  large 
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Figure  2.  Mean  velocity  profiles,  normalized  by  Uoo,  along  upper  surface.  Sym¬ 
bols:  case  A - ,  B - ,  measurements  by  Hastings  o  and  Wadcock  +  . 

streamwise  gradients  of  the  advected  variable.  This  situation  is  typically  encoun¬ 
tered  near  the  nose  and  the  trailing  edge  of  the  airfoil.  Other  sources  are  the  outflow 
boundary  (an  artificial  boundary  layer  is  generated  in  the  streamwise  direction)  and 
mesh  stretching.  As  shown  by  Cain  &  Bush  (1994),  waves  propagating  into  an  in¬ 
creasingly  coarse  (fine)  mesh  are  amplified  (dampened)  in  a  centered  scheme.  In 
our  simulation  we  find  that  strong  2-A-waves  appear  near  the  nose  and  near  the 
trailing  edge.  The  wiggles  appear  almost  exclusively  in  the  streamwise  coordinate 
direction.  Part  of  these  waves  travel  with  different  phase  speed  and  cancellation 
occurs.  However,  other  parts  are  steady  and  accumulate  in  time.  These  stand¬ 
ing  waves  contaminate  the  potential  flow  region  after  long  integration  times.  It  is 
difficult  to  assess  to  what  degree  the  solution  is  contaminated  by  the  presence  of 
2-A-waves.  On  a  staggered  mesh,  velocity  components  are  averaged  in  order  to 
obtain  fluxes  at  ceU  faces.  This  averaging  on  a  scale  of  the  mesh  cell  can  sometimes 
completely  hide  the  2- A- wave.  For  example,  the  convective  term  d(uu)/dx  in  the 
streamwise  momentum  balance  is  evaluated  as 
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The  finite  difference  expression  renders  the  same  value  independently  whether  an 
oscillatory  part  in  the  i-direction  Uj  =  (— !)*«„  with  zero  mean  and  arbitrary  am¬ 
plitude  Mo  is  added  or  not.  Similarly,  if  a  2-A-wave  in  the  i-direction  is  present 
in  the  v  velocity  component,  it  will  not  appear  in  the  discrete  approximation  for 
d{uv)/dy.  However,  it  will  contaminate  the  term  d{uv)ldx.  Time  averaged  fields 
of  velocity  components  show  2-A-waves  in  the  potential  flow  region,  but  the  pres¬ 
sure  field  is  virtually  free  of  wiggles.  This  indicates  that  the  presence  of  2-A-waves 
in  the  potential  flow  region  may  be  tolerated  to  a  certain  degree  since  wiggle  free 
streamlines  in  accordance  with  the  pressure  fleld  can  be  reconstructed. 
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Figure  3.  Pressiire  distribution  around  the  airfoil.  Symbols:  LES  ,  Wad- 

cock  o  . 

The  strongest  effect  of  2- A-waves  comes  from  the  associated  limitations  for  the 
computational  timestep.  Large  amplitude  wiggles  in  the  wall  normal  velocity  com¬ 
ponent  in  conjunction  with  rather  fee  wallnormal  spacing  cause  high  CFL  numbers 
near  the  nose.  The  resulting  timestep  limitations  are  so  severe  that  the  simulation 
can  not  be  carried  out  at  an  affordable  cost.  We  therefore  resorted  to  an  ad  hoc 
modification  of  the  numerical  scheme.  In  a  small  region  near  the  nose  (less  than 
2%  of  the  chord)  we  applied  a  l:2:l-filter  in  the  streamwise  and  spanwise  direction 
which  eflficiently  eliminates  all  2- A- waves.  Filtering  is  equivalent  to  adding  a  di¬ 
rection  dependent  diffusion  term  to  the  equations.  Justification  for  this  procedure 
comes  from  the  fact  that  the  flow  near  the  nose  is  laminar  and  filtering  on  a  scale 
of  the  grid  cell  does  not  affect  the  flow  physics.  Additionally,  the  boundary  layer 
in  the  experiments  was  tripped  at  a  location  around  x/c  =  0.02,  thereby  fixing  the 
region  of  laminar-turbulent  transition.  We  find  that  the  flow  spontaneously  transi¬ 
tions  as  soon  as  the  filter  ends.  In  this  sense,  we  control  the  location  of  transition 
by  setting  the  streamwise  extent  of  the  region  where  the  solution  is  filtered.  The 
filter  extended  about  40  layers  away  from  the  wall  and  faded  to  zero  over  another 
15  layers.  Unfortunately,  this  procedure  changes  the  potential  flow  sigmficantly. 
Because  the  mesh  cells  are  rather  large  in  the  outer  part  of  the  domain,  filtering  on 
the  grid  scale  is  no  longer  negligible  on  the  scale  where  the  potential  flow  changes 
near  the  nose.  Future  simulations  can  easily  avoid  this  problem  by  limiting  the  filter 
to  the  vicinity  of  the  surface,  i.e.  it  should  end  near  the  boundary  layer  edge.  No 
attempt  was  made  to  dampen  2-A-\^aves  in  the  trailing  edge  region  where  the  flow 
is  fully  turbtdent.  Any  filtering  there  would  probably  affect  the  flow  physics. 

2.4  Simulation  results  and  discussion 

Figure  1  shows  time  series  of  the  spanwise  velocity  fluctuation  w  recorded  at 
several  stations  along  the  upper  surface  of  the  airfoil.  We  observe  a  shift  in  the 
frequency  which  corresponds  to  the  most  energetic  motions  towards  lower  values  as 
the  recording  station  moves  closer  to  the  trailing  edge.  This  is  consistent  with  the 
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Figure  4.  Boundary  layer  thickness  8,  displacement  thickness  5*,  momentum 
thickness  9  and  shape  factor  H  =  8*  f  9  along  the  upper  surface  of  the  airfoil.  Sym¬ 
bols:  - LES,  •  Hastings,  x  Wadcock. 


increase  of  an  inertial  timescale  (ratio  of  the  boundary  layer  thickness  to  the  edge 
velocity)  as  the  boundary  layer  grows  under  the  influence  of  the  adverse  pressure 
gradient.  It  becomes  evident  that  the  solution  has  to  be  sampled  over  several  time 
units  c/Uoo  in  order  to  obtain  representative  turbulence  statistics  for  the  rear  part 
of  the  airfoil. 

Statistics  were  obtained  by  averaging  the  instantaneous  flow  fields  in  the  spanwise 
homogeneous  direction  and  in  time  over  more  than  2cfUoo-  Proflles  of  the  mean 
velocity  in  a  surface  normal  coordinate  system  are  shown  in  Fig.  2.  At  the  first  two 
stations,  the  edge  velocity  is  about  12%  smaller  than  measured  by  Hastings.  As 
mentioned  earlier,  this  is  a  side  effect  from  the  filter  which  was  applied  in  the  nose 
region  in  order  to  eliminate  2-A-waves.  Since  filtering  was  limited  to  a  region  close 
to  the  surface,  simulated  and  measured  mean  flow  agree  much  better  for  distances 
greater  than  yfc  =  0.06.  Although  a  better  match  between  simulated  and  measured 
edge  velocity  is  desirable  (and  can  easily  be  obtained  by  further  reducing  the  dis¬ 
tance  from  the  surface  over  which  the  filter  is  applied),  we  don’t  expect  turbulence 
statistics  to  be  significantly  affected.  One  reason  is  the  observation  that  the  simu¬ 
lated  adverse  pressure  gradient  matches  the  measured  one  over  most  of  the  upper 
surface,  see  Fig.  3.  Filtering  affects  mainly  the  magnitude  of  the  suction  peak  and  is 
partially  responsible  for  the  offset  in  the  simulated  pressure  distribution.  Addition¬ 
ally,  since  wind  tunnel  walls  were  not  properly  considered  in  this  simulation,  the 
pressure  distribution  near  the  nose  will  deviate  from  the  measured  one,  see  Jansen 
(1995).  The  goal  of  the  present  study  is  to  predict  the  boundary  layer  growth  and 
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Figure  5.  RMS  of  velocity  fluctuations  and  w',  normalized  with  Uoo- 

Symbols:  case  A - ,  B - ,  o  and  v'  □  from  Hastings.  In  terms  of  relative 

magnitude,  the  three  curves  for  each  simulation  are  v',  w',  and  «'  respectively. 

the  amount  of  separation  near  the  trailing  edge.  Accurate  prediction  of  the  suction 
peak  is  of  secondary  interest. 

Displacement  and  momentum  thickness  from  the  simulation  lie  in  between  the 
measurements  of  Hastings  &  Williams  (1984)  and  Wadcock  (1987)  upstream  of 
x/c  —  0.4,  see  Fig.  4.  The  experimental  v^ues  differ  by  up  to  40%  as  a  result  of  dif¬ 
ferences  in  boundary  layer  tripping  and  Reynolds  number.  However,  the  measured 
shape  factor  H  w  1.55  is  similar  in  both  experiments  in  the  region  r/c  =  0.2...0.4. 
Contrary  to  the  experiment,  H  drops  gradually  in  the  simulation  in  the  region 
xfc  =  0.2... 0.4  and  readies  values  as  low  as  1.4. 

Since  both  experiments  measure  similar  boundary  layer  growth  and  flow  retarda¬ 
tion  near  the  trailing  edge,  the  flow  development  does  not  seem  to  be  very  sensitive 
with  respect  to  the  exact  values  of  6*  and  6  of  the  turbulent  boundary  that  develops 
behind  the  transition  strip.  Although  the  thickness  of  the  simulated  boundary  layer 
is  close  to  the  measmred  ones  in  the  front  part  of  the  airfoil,  the  imderprediction 
of  the  shape  factor  in  the  simulation  and  the  initially  opposite  trend  (decline  as 
opposed  to  a  growth)  indicates  insufficiencies  in  the  simulated  boundary  layer  for  a 
considerable  part  of  the  upper  surface.  This  is  not  surprising  since  the  resolution  is 
so  coarse  that  the  near  wall  structures  can  hardly  be  resolved  properly.  Examina¬ 
tion  of  instantaneous  flow  fields  close  to  the  surface  reveals  a  very  streaky  structure 
with  typical  spacings  in  the  order  of  a  few  mesh  cells.  Similarly,  spanwise  two  point 
correlations  show  a  zero-crossing  within  2-3  spanwise  grid  points  for  all  near-wall 
locations  upstream  of  x/c  =  0.5,  see  Fig.  6.  This  indicates  that  the  simulation  has 
TnaryT^a.1  resolution  near  the  wall.  Further  evidence  comes  from  the  comparison  of 
the  present  case  with  an  earlier  simulation  which  was  a  factor  of  2  coarser  in  the 
streamwise  direction  and  a  factor  of  1.5  coarser  in  the  spanwise  direction.  The  flow 
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Figure  6.  Spanwise  two-point  correlations  Run  ( - )>  Rw  (  o  )  and  Ry,^ 

( -  )  versus  distance  2r/c  for  three  stations  along  the  upper  surface.  Bottom 

figures  correspond  to  the  near  wall  region,  top  figures  io  y  fS  m  0.3.  The  y  coordinate 
of  the  lower  figures  is  shifted,  i.e.  y  =  —1.5  corresponds  to  a  zero  correlation. 


retardation  and  the  boundary  layer  growth  was  significantly  improved  on  the  finer 
mesh.  Therefore,  further  grid  refinement  and,  subsequently,  a  better  prediction 
of  the  boimdary  layer  in  the  front  region  might  lead  to  better  agreement  between 
simulation  and  measurements  over  the  entire  airfoil. 

The  shear  stress  provided  by  the  SGS  model  is  an  indicator  for  the  role  of  the 
SGS  model.  The  maximum  contribution  is  about  15%  of  the  resolved  stress  xTo  and 
is  found  in  the  front  part  of  the  airfoil  where  the  resolution  is  coarse.  Near  the 
trailing  edge,  the  SGS  stress  is  negligible  compared  to  the  resolved  Reynolds  shear 
stress.  The  ratio  of  SGS  eddy  viscosity  to  molecular  viscosity  is  about  20,  which 
emphasizes  the  important  role  of  the  model  for  the  kinetic  energy  budget. 

RMS  values  of  the  velocity  fluctuations  are  shown  in  Fig.  5.  Agreement  between 
simulation  and  experiment  is  reasonable  in  the  middle  section  of  the  mrfoil.  In  a 
characteristic  manner  for  an  adverse  pressrure  gradient  boimdary  layer,  the  location 
of  maximiun  rms  values  (and  Reynolds  shear  stress)  moves  towards  the  outer  part  of 
the  boundary  layer.  Also,  the  anisotropy  of  the  fluctuations  in  the  outer  part  of  the 
boundary  layer  is  greatly  reduced.  Substantial  difference  between  simulation  and 
experiment  are  indicated  by  the  large  discrepancy  in  simulated  and  measured  rms 
values  (and  shear  stress)  near  the  trailing  edge.  It  is  imclear  whether  this  mismatch 
is  a  local  effect  or  rather  a  result  of  differences  in  (spatial)  flow  history  between 
experiment  and  simulation. 

Results  from  cases  A  and  B,  which  only  differ  with  respect  to  the  spanwise  domain 
size,  are  surprisingly  similar.  Two-point  correlations  from  the  outer  part  of  the 
boundary  layer  of  case  A  do  not  drop  to  zero  within  half  the  spanwise  width  for 
locations  downstream  of  x/c  =  0.6,  see  Fig.  6.  This  means  that  the  large  scales 
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of  motion  are  affected  by  the  presence  of  artificial  periodic  botmdaries.  Since  the 
limitations  are  much  more  severe  in  case  B  as  compared  to  A,  one  would  expect 
that  both  cases  deviate  in  the  rear  psirt.  Presently,  it  is  not  clear  why  the  simulation 
is  rather  insensitive  with  regard  to  the  domain  width.  Kaltenbach  (1994)  made  a 
similar  observation  for  a  flow  in  a  diflfuser  where  the  aspect  ratio  of  the  outlet  duct 
was  smaller  than  0.5.  Doubling  the  aspect  ratio  had  only  a  small  effect  on  the  flow 
evolution.  The  cost  for  case  B  is  about  half  that  of  case  A.  Further  studies  on  the 
effect  of  grid  refinement  would  be  much  cheaper  if  the  domain  width  of  case  B  turns 
out  to  be  sufficient. 

3.  Conclusions  and  future  goals 

Wall  resolving  LES  of  flow  aroxmd  an  airfoil  has  been  demonstrated  to  be  feasi¬ 
ble  with  present  computers  and  standard  numerical  schemes  for  LES.  Qualitatively, 
the  simulation  captures  typical  features  of  separating  flows  such  as  boundary  layer 
retardation  and  drastic  increase  in  Reynolds  stresses.  This  demonstrates  the  capa¬ 
bility  of  the  LES  concept  to  deal  with  flows  in  complex  configurations  of  immediate 
technical  interest.  However,  the  resolution  provided  was  probably  too  coarse  to 
adequately  simulate  the  boundary  layer  in  the  first  half  of  the  airfoil.  Although  the 
resolution  might  have  been  adequate  for  the  rear  part,  the  overall  agreement  with 
measurements  with  respect  to  prediction  of  backflow  is  not  satisfactory.  History 
effects  might  play  a  role,  and  further  studies  shotdd  attempt  to  match  better  the 
integral  boundary  layer  parameters  of  the  experiment  at  an  early  station.  Because 
of  conservation  properties,  the  use  of  centered  difference  schemes  is  very  desirable 
in  the  context  of  LES.  However,  the  emergence  of  2-A-waves  is  a  serious  problem 
for  the  present  high  Reynolds  number  flow  and  needs  further  consideration,  for  ex¬ 
ample,  usage  of  explicit  filters  as  explored  by  Lund  k  Kaltenbach  in  this  volume. 
Comparison  of  two  cases  with  different  domsin  width  did  not  show  significant  sen¬ 
sitivity  with  respect  to  this  parameter  in  the  range  considered.  Future  simulations 
should  consider  the  effect  of  wind  tuimel  (top  and  bottom)  walls  by  a  corresponding 
modification  of  domain  size  and  boimdary  conditions. 
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APPENDIX  3.  Conservative  properties  of 
finite  difference  schemes  for  incompressible  flow 

By  Youhei  MorinishP 


1.  Motivation  and  objectives 

The  ptirpose  of  this  research  is  to  construct  acciirate  finite  difference  schemes  for 
incompressible  unsteady  flow  simulations  such  as  LES  (large-eddy  simulation)  or 
DNS  (direct  numerical  simulation). 

Experience  has  shown  that  kinetic  energy  conservation  of  the  convective  terms 
is  required  for  stable  incompressible  tmsteady  flow  simulations.  Arakawa  (1966) 
showed  that  a  finite  difference  scheme  that  conserves  the  enstrophy  in  the  absence 
of  viscous  dissipation  is  required  for  long-time  integration  in  the  two-dimensional 
vorticity-streamfunction  formulation.  The  corresponding  conserved  variable  is  ki¬ 
netic  energy  in  velocity-pressure  formulation,  and  some  energy  conservative  finite 
difference  schemes  have  been  developed  for  the  Navier-Stokes  equations  in  three  di¬ 
mensions.  Staggered  grid  systems  are  usually  required  to  obtain  physically  correct 
pressure  fields.  The  standard  second  order  accurate  finite  difference  scheme  (Harlow 
&  Welch  1965)  in  a  staggered  grid  system  conserves  kinetic  energy  and  this  scheme 
has  proven  useful  for  LES  and  DNS.  However,  the  accuracy  of  the  second  order 
finite  difference  scheme  is  low  and  fine  meshes  are  required  (Ghosal  1995).  Spectral 
methods  (Canuto  et  al.  1988)  offer  supreme  accuracy,  but  these  methods  are  lim¬ 
ited  to  simple  flow  geometries.  Existing  fourth  order  accvurate  convective  schemes 
(A-Domis  1981,  Kajishima  1994)  for  staggered  grid  systems  do  not  conserve  kinetic 
energy.  Higher  order  staggered  grid  schemes  that  conserve  kinetic  energy  have  not 
been  presented  in  the  literature. 

The  conservation  of  kinetic  energy  is  a  consequence  of  the  Navier-Stokes  equations 
for  incompressible  flow  in  the  inviscid  limit.  In  contrast,  energy  conservation  in  a 
discrete  sense  is  not  a  consequence  of  momentum  and  mass  conservation.  It  is 
possible  to  derive  numerical  schemes  that  conserve  both  mass  and  momentum  but 
do  not  conserve  kinetic  energy.  It  is  also  possible  to  derive  schemes  that  conserve 
kinetic  energy  even  though  mass  or  momentum  conservation  axe  violated. 

In  this  report,  conservation  properties  of  the  continuity,  momentmn,  and  kinetic 
energy  equations  for  incompressible  flow  are  specified  as  analytical  requirements  for 
a  proper  set  of  discretized  equations.  Existing  finite  difference  schemes  in  staggered 
grid  systems  axe  checked  for  satisfaction  of  the  reqmrements.  Proper  higher  order 
accurate  finite  difference  schemes  in  a  staggered  grid  system  are  then  proposed. 
Plane  channel  flow  is  simulated  using  the  proposed  fourth  order  accurate  finite 
difference  scheme  and  the  results  compared  with  those  of  the  second  order  accurate 
Harlow  and  Welch  (1965)  algorithm. 

1  Permanent  address:  Nagoya  Institute  of  Technology,  Japan 
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2.  Accomplishments 

2.1  Analytical  requirements 

The  continiiity  and  momentum  equations  describe  the  motion  of  incompressible 
flow.  For  convenience  later  in  the  analysis,  these  equations  are  written  symbolically 
as 

{Cant.)  =  0  (1) 

d'o ' 

+  {Conv.)i  +  {Pres.)i  +  {Visc.)i  =  0  (2) 

where 

=  (pre».)i  =  ^,  (Visc.),  =  ^  (3).  (4),  (5) 

Here,  Vi  is  velocity  component,  p  is  pressure  divided  by  density,  and  Tij  is  viscous 
stress.  Henceforth,  p  will  be  referred  to  as  pressure. 

The  conservation  properties  of  Eqs.  (1)  and  (2)  will  now  be  estabhshed.  Note 
that  Eq.  (2)  is  in  the  following  form. 


•••  =0  (6) 

The  term  is  conservative  if  it  can  be  written  in  divergence  form 

=  V  •  (7) 

C/X  i 


To  see  that  the  divergence  form  is  conservative,  integrate  Eq.  (6)  over  the  volume 
and  make  use  of  Gauss’s  theorem  for  the  flux  terms  k  —  1, 2,  •  •  •,  all  of  which  are 
assumed  to  satisfy  Eq.  (7) 


(8) 


From  Eq.  (8),  we  notice  that  the  time  derivative  of  the  sum  of  ^  in  a  volume  V 
equals  the  sum  of  the  flux  on  the  surface  S  of  the  volume.  In  partictilar,  the 
sum  of  <j>  never  changes  in  periodic  field  if  is  conservative  for  all  k. 

Note  that  the  pressinre  (Pres.)j  and  viscous  terms  (V’*sc.)i  axe  conservative  a 
priori  in  the  momentmn  equation  since  they  appear  in  divergence  form.  The  con¬ 
vective  term  is  also  conservative  a  priori  if  it  is  cast  in  divergence  form.  This  is 
not  always  the  case,  however,  and  we  shall  investigate  alternative  formulations.  To 
perform  the  analysis,  we  regard  {Conv.)i  as  a  generic  form  of  the  convective  term 
in  the  momentum  equation.  At  least  four  types  of  convective  forms  have  been  used 
traditionally  in  analytical  or  numerical  studies.  These  forms  are  defined  as  follows. 

dvjVi 
dxj 


{Div.)i 


(9) 
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(10) 

_  1  dvjVi  1  1  _  dvi 
~  2  dxj  2^^  dxj 

(11) 

fdvi  5t;j\  \dvjVj 
dxi )  2  dxi 

(12) 

As  mentioned  above,  the  divergence  form,  (Div.)^  is  conservative  a  priori.  {Adv.)i, 
{Skew.)i,  and  (Rot.ji  are  referred  to  as  advective,  skew-symmetric,  and  rotational 
forms  respectively.  The  four  forms  are  connected  with  each  other  through  following 

{Adv.)i  =  {Div.)i  -  Vi  •  (Cont.)  (13) 


{Skew.)i  =  i(Diu.)i  +  ^{Adv.)i 
{Rot.)i  =  {Adv.)i 


(14) 

(15) 


We  notice  that  there  are  only  two  independent  convective  forms,  and  the  two  are 
equivalent  if  (Cont.)  =  0  is  satisfied.  It  is  also  apparent  that  the  advective,  skew- 
symmetric,  and  rotational  forms  are  conservative  as  long  as  the  continuity  equation 
is  satisfied. 

The  transport  equation  of  the  square  of  a  velocity  component,  /2,  is  Vj  times 
i  =  1  component  of  Eq.  (2). 


+  ^1 '  iConv.)\  +  vi  •  (Pres.)i  -1-  vi  •  (yisc.)i  =  0  (16) 

ot 

In  the  above  equation,  the  convective  term  can  be  modified  into  the  following  forms 
corresponding  to  those  in  the  momentum  equation. 


V, .  (Div.)i  =  ■  (.Coni.) 

(17) 

Vi  •  {Adv.)i  =  .  (Cont.) 

(18) 

.  dviVi^/2 

vi'{Skew.)i= 

(19) 

Note  that  the  skew-symmetric  form  is  conservative  a  priori  in  the  velocity  square 
equation.  Since  the  rotational  form  is  equivalent  to  advective  form,  the  four  con¬ 
vective  forms  are  conservative  if  (Cont.)  =  0  is  satisfied. 

The  terms  involving  pressure  and  viscous  stress  in  Eq.  (16)  can  be  modified  into 


following  forms. 


ui  •  (Pres.)i  = 


dpvi 

dxi 


(20) 
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Terms 

in  Momentum  Eq. 

TransportEquations 

vif2 

K 

(Div.) 

0 

0 

0 

{Adv.)  =  {Rot.) 

0 

0 

0 

{Skew.) 

0 

o 

0 

{Pres.) 

G 

X 

0 

(Vise.) 

O 

X 

X 

Table  1.  Conservative  properties  of  convective,  pressure,  and  viscous  terms  in 
the  Vi,  Ui/2,  and  K  equations.  0  is  conservative  a  'priori,  Q  is  conservative  if 
{Coni.)  =  0  is  satisfied,  and  x  is  not  conservative. 


vi  •  (yisc.)i 


drijVi  dvi 

dxj  dxj 


(21) 


These  terms  are  not  conservative  since  they  involve  the  pressure-strain  and  the 
viscous  dissipation. 

We  can  determine  the  conservative  properties  of  v^^  JT.  and  U3^/2  in  the  seune 
manner  as  for  12. 

The  transport  equation  of  kinetic  energy,  K  =  ViVi/2,  is  Vj  times  i  component  of 
Eq.  (2)  with  summation  over  i. 


dK 

+  Vi  •  {C(mv.)i  -H  Vi  •  (Pres.)i  +  Vi  ■  (Visc.)i  =  0 


(22) 


In  Eq.  (22),  the  conservation  property  of  the  convective  term  is  determined  in  the 
same  manner  as  for  V\^/2.  In  addition,  the  terms  involving  pressure  and  viscous 
stress  in  Eq.  (22)  can  be  modified  into  following  forms. 


Vi  •  {Pres.)i 


Vi  •  (Vise.),- 


dpvi 

dXi 

■  {Cont.) 

(23) 

dTijVi 

dxi 

dvi 
dx  i 

(24) 

The  pressure  term  in  Eq.  (22)  is  conservative  if  (Coni.)  =  0  is  satisfied.  The  viscous 
stress  term  in  Eq.  (22)  is  not  conservative  because  the  second  term  on  the  right-hand 
side  of  Eq.  (24)  is  the  energy  dissipation. 

Table  1  provides  a  summary  of  conservative  properties  of  convective,  pressure 
and  viscous  terms  in  the  transport  equations  of  Ui,  Vi/2  and  K  for  incompressible 
fiow.  The  final  goal  of  this  work  is  to  derive  higher  order  accurate  finite  difference 
schemes  that  satisfy  these  conservative  properties  in  a  discretized  sense. 


Discretized  operators 

Before  starting  the  main  discussion,  discretized  operators  need  to  be  defined.  In 
this  report,  the  discussion  of  the  discretized  equations  will  be  limited  to  uniform 
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grid  systems.  The  widths  of  the  numerical  grid  in  each  direction,  hi,  h2,  hi,  are 
constant.  The  grid  system  shown  in  Fig.  1  will  be  referred  to  as  a  staggered  grid 
system.  In  the  staggered  grid  system,  the  velocity  components  Ui  (i  =  1,2,3) 
are  distributed  around  the  pressure  points.  The  continuity  equation  is  discretized 
centered  at  pressure  points.  The  momentum  equation  corresponding  to  each  velocity 
component  is  centered  at  the  respective  velocity  point. 

Let  the  finite  difference  operator  acting  on  ^  with  respect  to  xi  and  with  stencil 
n  be  defined  as  follows. 

6nf>  _ ~Hy^hi/2,  X2,  3^3)  nhi/2,  X2,  3?$)  (25) 

SnXt  nhi 

xi,  *2,  xs 

Also,  an  interpolation  operator  acting  on  ^  in  the  xi  direction  with  stencil  n 

as  follows. 

^nxi  I  __  ^(^1  “b  nhi/2,  X2,  *^3)  "b  nhi  f^,  X2,  ^3)  (26) 

Ixi,  X2,  Xs  2 

In  addition,  define  a  special  interpolation  operator  of  the  product  between  <f>  and  if> 
in  the  xi  direction  with  stencil  n. 


=  -^(®i  +  n/n/2,  ®2>  xz)  V’(®1  ~  nhi/2,  X2,  xs) 

xi,  X2,  X8  2  ^27) 

+  +nhi/2,  X2,  xz)  <j){xi  —nhil2,  X2,  xz) 

A 

Equations  (25)  and  (26)  are  second  order  accurate  approximations  to  first  deriva¬ 
tive  and  interpolation,  respectively.  Combinations  of  the  discretized  operators  can 
be  used  to  make  higher  order  accurate  approximations  to  the  first  derivative  and 
interpolation.  For  example,  fomiih  order  accurate  approximations  are  as  follows. 

9  Si<f>  1  Sz<l>  ^  d<f>  3  <f>  4 

S6zXi~dxi  640  5a:i®  ^ 


(28) 
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9-jlXi  l-j3xi  ,  3  <j>  -  4  /nn\ 

Discretized  operators  in  the  X2  and  X3  directions  are  defined  in  the  same  way  as  for 
the  xi  direction. 

We  define  two  types  of  conservative  forms  in  the  discretized  systems.  in 

Eq.  (6)  is  (locally)  conservative  if  the  term  can  be  written  as 

(30) 

OlXj  02Xj  OzXj 

This  definition  corresponds  to  the  analytical  conservative  form  of  Eq.  (7).  is 
globally  conservative  if  the  following  relation  holds  in  a  periodic  field. 


*1  *2  *3 


The  sum  that  appears  in  Eq.  (31)  is  taJcen  over  the  period  of  respective  direction. 
AV  (=  hih2hz)  is  a  constant  in  a  uniform  grid  system.  The  definition  of  global 
conservation  corresponds  to  the  conservation  property  of  Eq.  (8)  in  a  periodic  field. 
The  condition  for  (local)  conservation  satisfies  the  condition  for  global  conservation. 

2.S  Continuity  and  pressure  term  in  a  staggered  grid  system 

Now  we  are  ready  to  consider  our  main  problem.  First  of  all,  let’s  examine  the 
conservative  property  of  the  pressure  term.  As  we  have  observed,  the  pressure  term 
should  be  conservative  in  the  transport  equations  of  momentum  and  kinetic  energy. 

In  the  staggered  grid  system,  define  the  discretized  continuity  and  pressure  term 
as  follows. 

(Cont.  —  S2)  =  =  0  (32) 

OiXi 

(Pres.  -  S2)i  =  ^  (33) 

diXi 

The  —52  denotes  that  the  above  approximations  axe  second  order  accurate  in  space. 
Fourth  order  approximations  for  the  continuity  and  pressure  term  in  the  staggered 
grid  system  are 

(Cant.  -  54)  0,  (34) 


(Pre«.-S4)iS§|l^-i|s^.  (33) 

^  '  SSiXi  BSzXi  ' 

Local  kinetic  energy  can  not  be  defined  uniquely  in  staggered  grid  systems  since  the 
velocity  components  zire  defined  on  staggered  grid  points.  Some  sort  of  interpolation 
must  be  used  in  order  to  obtain  the  three  components  of  the  kinetic  energy  at  the 
same  point.  The  required  interpolations  for  the  pressure  terms  in  the  and  K 
equations  are 


Sip  SiUip^^' 


■P'  (Cont  —  52), 
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• 

for  Momentum  Eq. 

Ui 

U?/2 

K 

{Pres.  -  52) 

0 

X 

Oi 

{Pres.  —  54) 

G 

X 

02 

Table  2.  Conservative  properties  of  finite  difference  schemes  for  the  pressure  term 
in  a  staggered  grid  system.  Q  is  conservative  a  priori,  0i  is  globally  conservative 
if  {Cant.  —  52)  =  0  is  satisfied,  O2  is  globally  conservative  if  {Cant.  —  54)  =  0  is 
satisfied,  and  x  is  not  conservative. 
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_  Irr  ^3P 

8  *  SiXi 


—  p  •  {Cant  —  54).  (37) 


8  8  S^xi  8  S^Xi  8  S^Xi 

The  following  relations  can  be  used  to  show  global  conservation  unambiguously. 


Xi  X2  ®3 


Xi  X2 


Xl  ®3  \  ^  /  Xl  X2  X3 

Therefore,  Eqs.  (33)  and  (35)  are  globally  conservative  if  the  corresponding  dis¬ 
cretized  continuity  equations  are  satisfied. 

Table  2  shows  the  summary  of  the  conservative  property  of  the  discretized  pres¬ 
sure  terms  in  a  staggered  grid  system. 

£.4  Convective  schemes  in  a  staggered  grid  system 

As  we  have  already  mentioned,  local  kinetic  energy  K  (=  UiUif2)  can  not  be 
defined  uniquely  in  a  staggered  grid  system.  Let  us  asstime  that  a  term  is  (locally) 
conservative  in  the  transport  equation  of  K  if  the  term  is  (locally)  conservative  in  the 
transport  equations  of  Ui^ /2,  t/2^/2  and  /2.  Since  the  conservative  properties  of 
172^12  and  Uz'^ / 2  are  estimated  in  the  same  manner  as  for  Ui^ /2,  only  conservative 
properties  of  convective  schemes  in  the  momentum  and  (2  equations  need  to  be 
considered. 

£.4.1  Proper  second  order  accurate  convective  schemes 

Define  second  order  accurate  convective  schemes  in  a  staggered  grid  system  as 

follows.  ,  , 

X  TT.*  Tr.» 

{Div.  -  52),  =  ^  \  —  (40) 


(Skew.  -  52)i  =  hoiv.  -  S2)i  +  -  52)j 


[ 
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FD  Schemes 
far  Momentum  Eq. 

TranspartEquations 

Ui 

Ut/2 

K 

{Div.  -  S2) 

0 

0 

0 

{Adv.  -  S2) 

0 

0 

0 

{Skew.  -  S2) 

0 

O 

G 

Table  3.  Conservative  properties  of  proper  second  order  accurate  convective  schemes  # 

in  a  staggered  grid  system.  Q)  is  conservative  a  priori  and  Q  conservative  if 
{Cant.  —  5'2)  =  0  is  satisfied. 


{Adv.  —  S2)i  is  connected  with  {Div.  —  5'2)i  through  the  following  relation. 

{Adv.  -  S2)i  =  {Div  -  S2)i  -  Ui  •  {Cont.  -  S2f^‘  (43) 


{Div.  —  S2)i  is  the  standard  divergence  form  in  a  staggered  grid  system  (Harlow  & 
Welch  1965).  {Adv.  —  S2)i  was  proposed  by  Kajishima  (1994).  {Skew.  —  S2)i  is 
equivalent  to  the  scheme  that  was  proposed  by  Pieicsek  &  Williams  (1970).  {Div.  — 
S2)i  is  conservative  a  priori  in  the  momentum  equation.  The  product  between  U\ 
and  {Skew.  —  S2)i  can  be  rewritten  as 

=  (44) 

OiXj 


Therefore,  {Skew.  —  S2)i  is  conservative  a  priori  in  the  transport  equation  of  Z7i^/2. 
By  using  Eq.  (43),  conservative  properties  of  the  various  schemes  are  determined. 
Table  3  shows  the  conservative  properties  of  {Div.—S2)i,  {Adv.—S2)i  and  {Skew.— 
S2)i.  These  schemes  are  seen  to  be  conservative  provided  continuity  is  satisfied.  In 
addition,  the  rotational  form  is  also  conservative  in  light  of  Eq.  (15). 


£..^.£  Proposal  of  proper  higher  order  accurate  convective  schemes 

It  is  of  interest  to  derive  a  proper  fourth  order  accurate  convective  scheme  for 
a  staggered  grid  system.  Existing  fourth  order  accurate  convective  schemes  for 
staggered  grid  systems  (  A-Domis  1981,  Kajishima  1994)  do  not  conserve  kinetic 
energy.  Here,  we  propose  the  following  set  of  fourth  order  accurate  convective 
schemes  in  a  staggered  grid  system. 


{Div.  —  S'4)i  = 


9  ^1 

8  SiXj 


1  ^3 

8  6zXj 


- 


(45) 


(46) 
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FD  Schemes 
for  Momentum  Eq. 

TransportEquations 

Ui 

Uf/2 

K 

{Div.  —  54) 

O 

0 

0 

{Adv.  —  54) 

0 

0 

o 

{Skew.  -  54) 

0 

o 

o 

l^ble  4.  Conservative  properties  of  proper  fourth  order  accurate  convective  schemes 
in  a  staggered  grid  system.  Q  is  conservative  a  priori  and  Q  is  conservative  if 
{Cant.  —  SA)  =  0  is  satisfied. 


(Skew.  —  S4)i  =  —{Div.  —  SA)i  +  —{Adv.  —  SA)i  (47) 

{Div.  —  SA)i  is  conservative  a  priori  in  the  momentum  equation.  The  product 
between  Ui  and  {Skew.  —  5'4)i  can  be  rewritten  as  follows. 


Ui  •  {Skew.  -  54)i  = 


-  W") 


1  ^3 

8  SzXj  \8  ^ 


Mi**' 
8^^  )  2 


(48) 


Thus,  {Skew.  —  5'4)i  is  conservative  a  priori  in  the  transport  equation  of  Ui^ /2. 
The  relation  between  {Adv.  -  54)<  and  {Div.  -  54),-  is  the  following. 


{Adv.  —  SA)i  —  {Div.  —  54)i  —  Ui  • 


-{Cant.-S4f^'  -  ^{Cant.  -  S4) 
8  8 


■3xi 


(49) 


This  equation  is  a  proper  discrete  analog  Eq.  (13),  and  {Adv.  —  S4)i,  {Div.  —  S4)i, 
and  {Skew.-S4)i  are  equivalent  if  (Coni. -54)  =  0  is  satisfied.  Using  this  relation, 
the  conservative  properties  of  the  present  schemes  are  determined.  Table  4  shows 
the  conservative  properties  of  the  present  schemes.  Comparing  Table  4  with  Table 
1,  we  see  that  the  present  schemes  are  a  proper  set  of  convective  schemes  provided 
that  the  continuity  equation  is  satisfied. 

Proper  higher  order  accurate  finite  difference  schemes  in  a  staggered  grid  system 
can  be  constructed  in  the  same  way  as  for  the  fourth  order  schemes. 

2.5  Channel  flow  simulation 

Numerical  tests  of  the  schemes  described  above  are  performed  using  plane  channel 
flow.  The  continuity  and  momentum  equations  for  incompressible  viscous  flow  are 
solved  using  the  proper  second  and  fourth  order  accurate  finite  difference  schemes 
in  a  staggered  grid  system  using  the  dynamic  subgrid  scale  model  (Germano  et  al. 
1991).  The  flow  is  drived  by  a  streamwise  pressure  gradient.  A  semi-implicit  time 
marching  algorithm  is  used  where  the  diffusion  terms  in  the  wall  normal  direction 


Youhei  Monnisht 


Figure  2.  LES  of  plane  channel  flow  at  Re=180  by  proper  second  and  fourth  order 
accurate  finite  difference,  (a)  Mean  streamwise  velocity;  (b)  Velocity  fluctuations. 

Symbols:  .  :  2nd  order  scheme; - :  4th  order  scheme;  •  :  DNS,  Kim,  et 

al  (1987); - ;  U+  =  5.5  =  2.5  log  y+. 
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are  treated  implicitly  with  the  Crank-Nicolson  scheme  and  a  third  order  Runge- 
Kutta  scheme  (Wray  1986)  is  used  for  all  other  terms.  The  fractional  step  method 
(Dukowicz  &  Dvinsky  1992)  is  used  in  conjimction  with  the  Van  Kan  (1986)  type 
of  pressure  term  and  wall  boimdary  treatment.  Periodic  boundary  conditions  are 
imposed  in  the  streamwise  and  spanwise  directions. 

The  subgrid-scale  model  is  the  d3mamic  model  (Germano  et  al.  1991)  with  the 
least  square  technique  (Lilly  1992).  Averaging  in  homogeneous  directions  is  used. 
Filtering  is  performed  in  the  spanwise  and  streamwise  directions. 

The  spatial  discretization  of  the  second  order  scheme  is  a  usual  one:  {Div.  —  S2) 
for  the  convective  term,  (Pres.  -  S2)  for  the  pressure  term,  and  (Cont.  -  S2)  for 
the  continuity.  The  corresponding  Poisson’s  equation  of  pressure  is  solved  using  a 
tri-diagonal  matrix  algorithm  in  wall  normal  direction  with  fast  Fourier  transforms 
(FFT)  in  the  periodic  directions.  The  second  order  accurate  control  volume  type 
^scretization  is  used  for  the  viscous  term. 

The  spacial  discretization  of  the  fourth  order  scheme  is  as  follows.  The  convec¬ 
tive  term,  the  pressure  term,  and  the  continuity  are  discretized  by  {Div.  —  54), 
(Pres.  -  54),  and  {Cont.  -  54),  respectively.  The  corresponding  Poisson’s  equa¬ 
tion  of  pressure  is  solved  using  a  septa-diagonal  matrix  algorithm  in  wall  normal 
direction  with  FFT  in  the  periodic  directions.  A  fourth  order  accurate  control  vol¬ 
ume  type  discretization  is  used  for  the  viscous  term.  The  subgrid  scale  terms  are 
estimated  with  second  order  finite  differences.  The  wall  boundary  condition  of  the 
fourth  order  scheme  is  designed  to  conserve  mass  and  momentum  in  the  wall  normal 
direction  in  a  discretized  sense. 

The  Reynolds  number  based  on  channel  half  width  and  wall  friction  velocity.  Re, 
is  180.  The  computational  box  is  47r  x  2  x  |7r,  and  the  mesh  contains  32  x  65  x  32 
points  (streamwise,  wall-normal,  and  spanwise  respectively). 

Figure  2  shows  the  profiles  of  mean  streamwise  velocity  and  velocity  fiuctuations 
from  the  proper  second  and  fomth  order  schemes.  Filtered  DNS  data  (Kim  et  al. 
1987)  are  plotted  as  a  reference  in  the  figures.  The  mean  streamwise  velocity  profile 
from  the  second  order  scheme  is  shifted  up  in  the  logarithmic  region.  This  defect  of 
the  second  order  scheme  is  usually  observed  in  coarse  LES  (Cabot  1994).  Another 
defect  of  the  second  order  scheme  in  coarse  LES  is  the  peak  value  of  streamwise 
velocity  fluctuation  is  too  high  (Cabot  1994).  These  defects  are  improved  by  using 
the  fourth  order  scheme.  The  computational  cost  of  the  fourth  order  method  is 
about  1.9  times  that  for  the  second  order  method. 

3.  Future  plans 

The  fourth  order  scheme  will  be  tested  in  high  Reynolds  number  channel  flow 
to  see  if  it  has  a  greater  advantage  when  the  velocity  fluctuations  have  a  relatively 
larger  fraction  of  energy  near  the  cutoff  waventunber. 
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APPENDIX  4.  Large-eddy  simulation 
of  flow  around  a  NACA  4412 
airfoil  using  unstructured  grids 

By  Kenneth  Jansen^ 


1.  Motivation  and  objectives 

Large-eddy  simulation  (LES)  has  matured  to  the  point  where  application  to  com¬ 
plex  flows  is  desirable.  The  extension  to  higher  Rejmolds  numbers  leads  to  an  im¬ 
practical  number  of  grid  points  with  existing  structured-grid  methods.  Furthermore, 
most  real  world  flows  are  rather  difficult  to  represent  geometrically  with  structured 
grids.  Unstructured-grid  methods  offer  a  release  from  both  of  these  constraints. 
However,  just  as  it  took  many  years  for  structured-grid  methods  to  be  well  im- 
derstood  and  reliable  tools  for  LES,  imstructured-grid  methods  must  be  carefully 
studied  before  we  can  expect  them  to  attain  their  full  potential. 

In  the  past  three  years,  important  building  blocks  have  been  put  into  place,  mak¬ 
ing  possible  a  careful  study  of  LES  on  unstructured  grids.  The  first  building  block 
was  an  efficient  mesh  generator  which  allowed  the  placement  of  points  according  to 
smooth  variation  of  physical  length  scales.  This  variation  of  length  scales  is  in  all 
three  directions  independently,  which  allows  a  large  reduction  in  points  when  com¬ 
pared  to  structured-grid  methods,  which  can  only  vary  length  scales  in  one  direction 
at  a  time.  The  second  building  block  was  the  development  of  a  dynamic  model  ap¬ 
propriate  for  unstructured  grids.  The  principle  obstacle  was  the  development  of 
an  unstructured-grid  filtering  operator.  New  filtering  operators  were  developed  in 
Jansen  (1994).  In  the  past  year,  some  of  these  filters  have  been  implemented  into 
a  highly  parallelized  finite  element  code  based  on  the  Galerkin/least-squares  finite 
element  method  (see  Jansen  et  al.  1993  and  Johan  et  al.  1992). 

We  have  chosen  the  NACA  4412  airfoil  at  maximtun  lift  as  the  first  simulation 
for  a  variety  of  reasons.  First,  it  is  a  problem  of  significant  interest  since  it  would 
be  the  first  LES  of  an  aircraft  component.  Second,  this  flow  has  been  the  subject  of 
three  experimental  studies  (Coles  and  Wadcock  1979,  Hasting  and  W^illiams  1987, 
and  Wadcock  1987).  The  ffist  study  found  the  maximtun  lift  angle  to  be  13.87®. 
The  later  studies  found  the  angle  to  be  12®.  Wadcock  reports  in  the  later  study 
that  the  early  data  agree  very  well  with  his  new  data  at  12®,  suggesting  that  the 
early  experiment  suffered  from  a  non-parallel  mean  flow  in  the  Caltech  wind  timnel. 
It  should  be  pointed  out  that  the  Reynolds-averaged  simulations  are  usually  nm 
at  13.87®  and  do  not  agree  with  the  data  when  r\m  at  12®  as  was  shown  in  Jansen 
(1995).  It  is  hoped  that  LES  can  clarify  this  controversy.  The  third  reason  for 
considering  this  flow  is  the  variety  of  flow  featiues  which  provide  an  important  test 
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Figure  1.  A  transition  strip  is  modeled  geometrically  by  applying  a  no-slip 
boundary  condition  to  the  nodes  which  form  a  surface  of  height,  shape,  and  position 
equivalent  to  Wadcock’s  serrated  tape  which  was  applied  to  the  airfoil  surface. 


of  the  dynamic  model.  Starting  from  the  nose  where  the  flow  stagnates,  thin  laminar 
boundary  layers  are  formed  in  a  very  favorable  pressure  gradient.  This  pressture 
gradient  soon  turns  adverse,  driving  the  flow  toward  a  leading  edge  separation. 
Only  the  onset  of  turbulence  can  cause  the  flow  to  remain  attached  or  to  reattach 
if  it  did  separate.  The  persistent  adverse  pressure  gradient  eventually  drives  the 
turbulent  flow  to  separate  in  the  last  20  percent  of  chord.  The  separation  bubble  is 
closed  near  the  trailing  edge  as  the  retarded  upper  surface  boundary  layer  interacts 
with  the  very  thin  lower  surface  boxmdary  layer.  The  large  difference  in  botmdary 
layers  creates  a  challenging  wake  to  simulate.  Only  the  dynamic  model  can  be 
expected  to  perform  satisfactorily  in  this  variety  of  situations;  from  the  laminar 
regions  where  it  must  not  modify  the  flow  at  all  to  the  turbulent  boimdary  layers 
and  wake  where  it  must  represent  a  wide  variety  of  subgrid-scale  structures. 

The  flow  configuration  we  have  chosen  is  that  of  Wadcock  (1987)  at  Reynolds 
number  based  on  chord  Rec  —  Uqoc/v  =  1.64  x  10®,  Mach  number  M  =  0.2,  and 
12**  angle  of  attack. 
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2.  Accomplishments 

2.1  Effect  of  wind  tunnel  walls  and  transition  strip 

In  Jansen  (1995)  a  grid  independent  solution  was  obtained  which  did  not  agree 
well  with  the  experiments  in  the  separated  region.  This  was  not  completely  sur¬ 
prising  since  two  important  effects  of  the  experiment  were  not  accounted  for  in  the 
simulation:  the  wind  tunnel  walls  and  the  transition  strip. 

Wadcock  used  a  strip  of  tape  with  serrations  cut  into  the  edge  on  the  upstream 
side.  The  serrated  tape  has  been  modeled  in  a  coarse  fashion  by  om  current  simu¬ 
lation  as  can  be  seen  in  Fig.  1.  The  tape  is  effectively  a  forward  facing  step  (with 
serrations)  of  height  ^99  / 4,  followed  by  a  baricward  facing  step. 

The  blockage  effect  of  the  wind  tunnel  walls  has  also  been  included  in  the  recent 
calculations.  Note  that  the  bovmdary  layers  on  the  walls  are  not  simulated;  rather, 
slip  boimdary  conditions  are  applied  on  the  wind  tuimel  walls  as  can  be  seen  in 
Fig.  2. 

These  two  effects  were  studied  separately  for  a  short  period  of  time  (not  sufl&cient 
for  converged  statistics  in  the  trailing  edge  region)  and  agreement  with  experiments 
was  seen  to  improve  in  both  cases.  The  effect  of  the  walls  was  somewhat  greater 
than  that  of  the  transition  strip.  This  discussion  is  left  qualitative  because  the 
enormous  cost  of  these  calculations  led  us  to  abandon  the  individual  effect  studies 
in  favor  of  using  our  limited  resources  on  converging  the  combined  effect  simulation. 
The  results  of  this  simulation  can  be  seen  in  Fig.  3  where  the  velocity  profiles  of 
the  new  simulation  (with  wind  timnel  walls  and  transition  strip)  are  compared  to 
the  original  simulation.  Note  the  large  increase  in  the  degree  of  separation. 

Though  the  new  simulation  is  in  better  agreement  with  the  experiment,  new 
problems  were  created.  As  can  be  seen  in  Fig.  4,  there  is  a  reduction  in  the  three- 
dimensionality  of  the  tmbulence  in  the  separated  region.  A  similar  story  is  told 
by  the  two-point  correlation  in  this  region,  which  shows  very  little  decay.  While 
some  reduction  is  to  be  expected,  concern  developed  as  to  whether  the  periodic 
boundary  condition,  which  is  applied  in  the  spanwise  direction,  was  promoting 
spanwise  coherent  vortices  due  to  insufficient  spanwise  extent.  Prom  Pig.  3  it  is 
apparent  that  the  boimdary  layer  is  significantly  thicker  in  the  new  simulation. 
With  a  spanwise  domain  width,  W  of  2.5%  of  chord  the  spanwise  domain  becomes 
less  than  a  boxmdary  layer  height  at  about  two-thirds  of  the  chord  length. 

2.2  Wider  domain  simulation 

The  most  obvious  choice  of  doubling  the  domain  while  maintaining  the  current 
resolution  was  postponed  in  lieu  of  a  doubling  of  the  domain  by  doubling  the  span- 
wise  size  of  each  element.  The  rational  for  this  decision  was  that  the  refinement 
studies  of  Jansen  (1995)  showed  only  a  small  change  from  the  current  grid  to  the 
twice  coarser  grid.  The  new  simulation  would  also  not  engender  an  increased  cost 
since  the  number  of  nodes  remained  the  same.  It  was  assumed  that  changes  in 
the  first  half  of  the  airfoil  would  suggest  inadequate  resolution  while  changes  in  the 
trailing  edge  region  would  address  the  domain  width  question. 
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Figure  2.  The  cross  sectional  plane  of  an  unstructured  mesh  which  accounts  for 
the  inviscid  effects  of  the  wind  tunnel  walls. 


tt 
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Figure  3.  Profiles  of  tangential  velocity  component  at  various  positions  along  the 
airfoil  surface  (r/c  =  0.59,0.66,0.78,0.82,0.95).  Solutions  correspond  to:  without 

wind  tunnel  walls  or  transition  strip - —  ,  with  wind  timnel  walls  and  transition 

strip - ,  Wadcock  □  ,  Hastings  and  Williams  o  . 


The  times-series  firom  the  new  simulation  are  presented  in  Fig.  5  where  the  three- 
dimensionality  can  be  seen  to  return  to  the  separated  region.  The  amount  of  sep¬ 
aration  is  reduced,  causing  some  departure  from  the  experimental  data  as  can  be 
seen  in  Fig.  6.  This  is  to  be  expected  because  the  presence  of  three-dimensional 
vortical  structures  in  the  separated  region  pump  high  momentum  fluid  down  to  the 
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Figure  4.  Time  series  of  the  spanwise  velocity  fluctuation  at  various  positions 
along  the  chord  length,  approximately  half  a  boundary  layer  height  off  of  the  wall. 
From  the  bottom  up  (xfc  =  0.1,0.3,0.529,0.66,0.815,0.95).  Note  that  the  top 
ctuwe  (ar/c  =  0.95)  indicates  a  loss  of  three-dimensionality  in  the  separated  region. 

wall,  reducing  the  magnitude  of  the  separation.  Fig.  6  also  illustrates  that  the  flrst 
half  of  the  airfoil  is  nearly  grid  independent. 

3.  Future  plans 

3.1  More  spanwise  domain  studies 

Based  on  the  fiuHings  in  the  previous  section  more  attention  will  be  given  to 
the  spanwise  domain  effects.  Since  the  resolution  changed  at  the  same  time  as  the 
expansion  of  the  domain,  it  is  difficult  to  isolate  the  two  effects.  For  this  reason  the 
first  study  will  maintain  the  5%  chord  domain  width  and  improve  the  resolution  in 
the  second  half  of  the  airfoil  where  the  solution  has  shown  some  change.  A  grid  has 
already  been  generated  to  accomplish  this  task,  and  a  simulation  has  just  begun. 
The  number  of  nodes  has  gone  up  by  70%,  making  this  simulation  significantly  more 
expensive  than  the  others  described  in  this  report. 
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Figure  5.  Time  series  of  the  spanwise  velocity  fluctuation  as  described  in  Fig.  4. 
Note  that  the  top  curve  (x/c  =  0.95)  has  a  strong  signal  indicating  a  return  of 
three-dimensionality  in  the  separated  region  with  the  increased  domain  width. 

S.2  Higher  order  methods 

Given  the  number  of  points  that  are  required  to  obtain  a  grid-independent  solu¬ 
tion,  it  seems  clear  that  higher  order  methods  should  be  explored.  This  is  straight¬ 
forward,  but  non-trivial,  to  do  with  the  finite  element  method.  There  are  two 
benefits  to  higher  order  methods  besides  the  obvious  one  of  higher  accuracy.  First, 
the  higher  order  methods  will  have  a  more  complete  representation  of  the  residual 
error  of  the  discrete  approximation  and,  therefore,  the  scheme  will  be  less  dissipa¬ 
tive.  Second,  alternative  filters,  described  in  Jansen  (1994),  can  be  implemented  and 
tested.  It  is  difficult  to  predict  at  this  time  if  the  method  will  lose  computational 
eflBciency  when  extended  to  higher  order. 

3.3  Computational  platform  change 

In  the  past  year  the  code  has  been  ported  to  the  IBM  SP2  (see  Bastin  1996  for 
details).  This  port  involved  the  use  of  MPI,  a  communication  standard  that  is  more 
widely  used  than  that  of  the  original  code,  which  should  make  the  port  to  other 
platforms  reasonably  simple.  In  the  coming  year  more  effort  will  be  appfied  in  this 
area  to  try  to  take  advantage  of  the  changing  computational  resources  available  for 
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Figure  6.  Profiles  of  tangential  velocity  component  at  various  positions  along  the 
airfoil  surface  {x/c  =  0.59,0.066,0.78,0.82,0.95).  Solutions  correspond  to:  without 

wind  ttumel  walls  or  transition  strip - ,  with  wind  tunnel  walls  and  tr^sition 

strip  ( W/c  =  0.025) - ,  with  wind  tunnel  walls  and  transition  strip  {W/c  =  0.05) 

- ,  Wadcock  □  ,  Hastings  and  Williams  o  . 

this  type  of  simulation,  thereby  expediting  its  progress. 
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